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Abstract 

In  this  paper  we  will  present  early  work  on  using  a  DDDAS  based  approach  to  the  construction  of  probabilistic 
estimates  of  volcanic  ash  transport  and  dispersal.  Our  primary  modeling  tools  will  be  a  combination  of  a  plume 
eruption  model  BENT  and  the  ash  transport  model  PUFF.  Data  from  satellite  imagery,  observation  of  vent  parameters 
and  windfields  will  drive  our  simulations.  We  will  use  uncertainty  quantification  methodology  -  polynomial  chaos 
quadrature  in  combination  with  data  integration  to  complete  the  DDDAS  loop. 
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1.  Introduction 

The  dynamic  data  driven  application  system  paradigm  (DDDAS)  is  formally  defined  in  a  series  of  recent  reports 

[1,  2]  and  in  earlier  reports  on  the  subject  as  systems  where  data  are  dynamically  integrated  into  simulation(s)  to 

augment  or  complement  the  application  model,  and,  where  conversely  the  executing  simulation  steers  the  measure¬ 
ment  (instrumentation  and  control)  processes  of  the  application  system.  The  full  range  of  research  issues  raised  by 
the  DDDAS  paradigm  are  many  and  diverse.  We  describe  here  the  construction  of  a  framework  for  DDDAS  to  target 
predictions  of  ash  clouds  that  are  produced  by  explosive  eruptions  from  volcanoes  and  address  primarily  questions 
related  to  the  quantification  of  uncertainty  and  a  computational  strategy  for  data  driven  forecasts  using  large  scale 
computing  and  analytics.  The  motivation  of  the  research  proposed  here  will  be  to  develop  techniques  that  can  be 
completed  in  a  timely  fashion  to  enable  the  DDDAS  control  loop  on  the  simulation  and  observation,  thus,  greatly 
enhancing  the  quality  of  the  prediction.  Recent  eruptions  of  Eyjafjallajokullin  Iceland  and  Puyehue  in  Chile  illustrate 
the  potential  of  these  events  to  interrupt  aviation.  Our  goal  is  to  provide  timely  and  accurate  ash  cloud  forecasts  to 
minimize  these  disruptions. 
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Our  core  approach  is  to  use  weighted  ensembles  of  a  novel  coupling  of  an  eruption  column  model  bent  and  the 
popular  puff  tool  as  our  simulation  tool  inside  a  DDDAS  loop.  Infrared  and  visual  satellite  imagery  can  provide  data 
on  the  position  of  these  clouds.  However,  both  simulation  and  observation  are  subject  to  large  error  and  uncertainty. 
Thus,  we  will  develop  methodology  to  carefully  characterize  the  error  and  uncertainties  in  a  timely  fashion.  Our 
proposed  techniques,  if  successful  will  greatly  improve  the  quality  of  information  available  on  the  potential  location 
and  concentrations  of  the  ash  clouds.  In  this  paper,  we  provide  an  overview  of  the  current  and  planned  efforts  and 
preliminary  results  in  implementing  a  DDDAS  approach  to  this  problem  domain.  In  a  companion  paper  [3]  we  outline 
in  detail  one  of  the  critical  phases  -  fast  solution  of  the  inverse  problem  to  identify  and  update  the  inputs  based  on 
observations. 

2.  Background 

2.7.  Volcanic  Plumes 

Long-range  volcanic  ash  transport  models  [4]  have  been  given  the  general  name  Volcanic  Ash  Transport  and 
Dispersion  Models  (VATDMs).  Numerous  VATDMs  have  been  developed  by  both  civil  and  military  aviation  or  me¬ 
teorological  agencies  to  provide  forecasts  of  ash  cloud  motion  [5].  The  use  of  these  models  became  well-known 
following  the  eruption  of  Eyjafjallajokull,  Iceland,  in  April,  2010,  when  the  NAME  model  [6]  forecasts  were  particu¬ 
larly  prominent  in  being  used  as  the  basis  for  an  almost  complete  shutdown  of  European  air  traffic.  In  this  paper,  we 
consider  a  VATDM  named  puf  f  that  is  used  by  the  Anchorage  Volcanic  Ash  Advisory  Center  ( VAAC)  for  real-time 
civil  aviation  forecasts.  Tanaka  [7]  and  Searcy  et  al.  [8]  developed  puff  for  predicting  the  paths  of  volcanic  clouds 
at  meso-  and  synoptic  scales,  puff  simplifies  the  eruption  plume  to  a  vertical  source,  and  uses  a  Lagrangian  pseudo¬ 
particle  representation  of  the  ash  cloud  in  a  detailed  3-D  regional  windfield  to  determine  the  trajectory  of  the  cloud, 
puf  f  and  other  dispersion  models  have  proven  extremely  useful  in  modeling  the  distal  transport  of  ash  for  aviation 
safety  [8,  5].  While  the  work  reported  here  is  in  the  context  of  puff  the  statistical  and  stochastic  methodology  that 
will  be  developed  is  independent  of  the  puff  model,  and  could  be  directly  applied  using  other  dispersion  models.  We 
are  in  the  process  of  implementing  this  methodology  with  the  more  recently  formulated  WRF-CHEM  tool. 

2.2.  puf  f  Simulation  Model 

There  are  over  a  dozen  operational  VATDMs  in  existence  that  have  been  used  for  ash  cloud  forecasting  and  post¬ 
event  analysis  by  the  world’s  nine  VAACs  [9,  5,  10,  11].  puff  is  one  of  these  VATDMs,  and  is  particularly  tied  to 
usage  by  the  Anchorage  VAAC.  The  model  is  intended  to  be  a  rapid  response  prediction  tool  of  modest  complexity, 
that  runs  on  workstation-class  computers.  Unlike  other  VATD  models,  it  is  designed  solely  for  predictions  of  ash-cloud 
movement  in  space  and  time,  and  it  generates  displays  that  can  be  tailored  to  user  needs  at  volcano  observatories,  in 
a  real-time  setting  during  an  eruption  crisis.  For  these  reasons, puf  f  is  necessarily  simple  by  design,  but  it  can  be 
intricate  in  its  implementation.  The  model  takes  into  account  the  predominant  physical  processes  that  control  particle 
movement,  such  as  winds,  dispersion  and  settling,  but  it  does  not  account  for  small  scale  physical  processes  that  play 
a  lesser  role.  In  some  cases  we  do  not  understand  in  detail  how  these  sub-scale  process  impact  the  overall  cloud 
movement  (e.g.  inter-particle  interactions),  in  other  cases  the  required  information  may  not  be  available  in  a  real-time 
situation  (e.g.  precipitation),  and  in  yet  other  cases  inclusion  of  the  process  may  significantly  lengthen  model  run-time 
(e.g.  radiative  heating). 

During  an  eruption  crisis,  puff  predictions  have  been  used  to  estimate  ash  cloud  movement  critical  to  the  assess¬ 
ment  of  potential  impacts  -  for  example,  on  aircraft  flight  paths,  or  on  ash  fall  at  airports,  or  for  health  concerns.  As 
other  more  complicated  (and  thus  slower)  VATD  model  results  become  available  they  can  be  incorporated  into  an 
impact  analysis  by  volcanologists  and  atmospheric  scientists.  Thus,  uncertainty  analysis  and  data  driven  correction  of 
VATD  models  in  general,  and  of  the  puf  f  model  in  particular,  is  important  for  appropriate  response  to  the  resulting 
hazards. 

puff  can  be  run  using  one  of  several  numerical  weather  prediction  (NWP)  windfields  [12-15].  These  NWP 
models  are  available  at  differing  levels  of  spatial  and  temporal  resolution;  local  use  of  the  WRF  (Weather  Research 
Forecast)  model  will  allow  for  customizing  the  model  domains  used  for  verification  and  validation,  puff  tracks  a 
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finite  number  of  Lagrangian  point  particles  of  different  sizes,  whose  location  R  is  propagated  from  timestep  k  to 
timestep  k  +  1  via  an  advection/diffusion  equation 

Ri(tk+ 1)  -  Ri(tk)  +  W(tk)At  +  Z(tk)At  +  S i(tk)At  (1) 

Here  R*(4)  is  the  position  vector  of  the  ith  particle  at  time  kAt,  W(tk )  is  the  local  wind  velocity  at  the  location  of  the 
ith  particle,  Z(4)  is  a  turbulent  diffusion  that  is  modeled  as  a  random  walk,  and  Sj(4)  is  a  source  term  which  models 
the  fallout  of  the  Ith  particle  due  to  gravity. 

Ash  loading  and  plume  height  are  thought  to  be  the  most  important  input  variables  for  volcanic  ash  transport  and 
dispersion  (VATD)  modeling,  such  as  that  applied  to  the  Eyjafjallajokull  plume  [16].  Current  usage  is  to  derive  mass 
loading  from  the  product  of  eruption  duration  and  mass  eruption  rate  (MER)  (equivalent  to  mass  flux,  Q ).  MER  in 
turn  is  typically  calculated  from  an  empirical  relationship  derived  for  strong  plumes  (w  »  v),  e.g.,  HT  =  1.67g0,259 
[17],  where  w  is  characteristic  plume  speed  (m/s),  v  is  wind  speed,  HT  is  eruption  column  height  (km)  and  Q  is 
mass  eruption  rate  given  as  the  equivalent  volume  eruption  rate  (cu  m/s)  assuming  a  magmatic  density  (for  e.g.  3000 
kg/cu  m  for  basalt  in  the  case  of  Eyjafjallajokull).  This  yields  ash  loading  as  a  function  of  cloud  height,  HT  alone. 
Plume  height,  however,  is  a  complex  function  of  source  and  environmental  conditions,  so  using  plume  height  alone  to 
estimate  mass  loading  can  dramatically  underestimate  mass  loading  in  the  case  of  high  wind. 

2.3.  bent  Eruption  Column  Model 

To  deal  with  overdependence  on  plume  height,  we  can  use  measurements  of  several  more  primitive  source  vari¬ 
ables:  vent  speeds  and  vent  radii,  and  ground  measurements  of  ash  properties.  Using  these  primitive  variables  as 
inputs  to  suitable  models  we  postulate  that  we  can  better  sample  the  uncertain  input  parameter  space  for  a  VATD 
model  rather  than  the  current  practice  of  using  plume  height  conversion. 

We  thus  track  ash  dispersal  using  an  eruption  column,  plume  rise  model  bent  with  uncertain  input  parameters 
of  vent  speed  and  radius,  and  mean  and  variance  of  grain  size,  coupled  with  a  radiosonde  from  near  the  initiation  of 
the  eruption  when  available  to  construct  a  sample  space  of  uncertain  input  parameters  of  grain  size,  mass  loading  and 
injection  height.  The  mass  loading  and  height  from  the  plume  rise  model  are  then  used  as  input  to  a  VATD  model, 
using  NCEP  Reanalysis  windfields  to  track  the  propagation  of  ash. 

The  bent  integral  eruption  column  model  will  be  used  to  produce  eruption  column  parameters  (mass  loading, 
column  height,  grain  size  distribution)  given  a  specific  atmospheric  sounding  and  source  conditions  [18].  bent  takes 
into  consideration  atmospheric  (wind)  conditions  as  given  by  atmospheric  sounding  data.  Thus  plume  rise  height  is 
given  as  a  function  of  volcanic  source  and  environmental  conditions. 

3.  DDDAS  Workflow  Characterization 

We  begin  by  laying  out  the  overall  workflow  and  use  it  to  describe  the  DDDAS  framework  for  this  application. 
Figure  1  describes  the  overall  workflow  of  the  ash  plume  prediction  DDDAS  application.  Key  to  workflow  are  the 
propagation  of  uncertainty  -  in  parameters  and  source  terms  (wind  field)  and  integration  with  data  from  satellites  and 
other  observations.  We  start  with  a  description  of  the  pdf  of  the  concentration  of  ash  at  initial  time  at  all  points  on  some 
grid  covering  the  jurisdiction  of  interest .  As  time  evolves  the  concentration  pdf  is  propagated  by  using  the  uncertainty 
described  by  the  wind  and  input  parameters.  The  wind  effect  is  accounted  for  by  using  AGMM  (See  details  below) 
in  the  inner  loop  in  lock-step  with  the  bent-puff  coupled  model.  The  outer  PCQ  loop  (see  Fig.  1)  accounts  for 
the  input  parameter  uncertainty.  The  propagated  concentration  pdf  is  updated  using  data  assimilation  methodology 
as  described  below.  After  every  data  assimilation  step  the  input  model  parameters  pdf  is  updated  using  the  Bayesian 
framework.  Based  on  the  modified  pdf  a  new  PCQ  based  sampling  must  be  redone.  Note  that  we  can  always  start  the 
process  at  some  posterior  time  if  a  concentration  field  and  its  variability  are  available.  Depending  on  the  uncertainty 
associated  with  the  posterior  pdfs  we  can  use  higher  fidelity  models  in  the  data  assimilation  loop. 

3.1.  Uncertainty  and  Stochastic  Variability 

There  are  a  number  of  significant  challenges  associated  with  uncertainty  characterization  and  propagation  for 
spatio-temporal  varying  systems  like  those  encountered  in  the  ash  plume  transport.  Uncertainty  arises  from  several 
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BENT:  Eruption  Plume  Model 
PUFF:  Ash  transport  and 
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sensors  for  ash  detection 


Figure  1:  Planned  workflow  of  this  DDDAS  application  as  mapped  to  a  set  of  heterogeneous  resources.  Blue  arrows  denote 
components  that  have  been  implemented,  while  brown  arrows  denote  work  in  progress. 


sources:  parameters,  initial  and  boundary  conditions,  forcing  functions,  are  known  only  to  certain  precision.  For 
example,  bent -puff  input  parameters  include  coefficient  of  gravity  fallout,  vent  diameter,  vent  velocity,  mean  grain 
size  and  windfield.  All  of  these  parameters  are  known  only  approximately  and  only  to  a  fixed  fidelity.  Despite  the 
potential  risk  to  property  and  life  from  inaccurate  prediction  of  ash  clouds,  there  has  never  been  a  thorough  quantita¬ 
tive  assessment  of  ash  cloud  predictions  due  to  parametric  and  stochastic  uncertainties.  Most  of  the  literature  [19-23] 
focuses  on  tracking  the  ash  material  plume  while  making  use  of  the  sensitivity  information  of  plume  location  to  mea¬ 
sured  material  flux  or  concentration.  While  a  detailed  sensitivity  analysis  can  relate  the  variations  in  input  parameters 
to  ash  cloud,  uncertainty  analysis  casts  a  much  broader  net  in  terms  of  assessing  confidence  of  predictions  based  on 
all  available  information.  These  aforementioned  two  questions  clarify  the  main  outcome  of  our  work  here  To  provide 
a  prediction  of  ash  cloud  motion,  together  with  a  quantitative  measure  of  the  variability  of  that  prediction,  informed  by 
all  available  data  in  real  time.  [Note:  By  “real  time”,  we  mean  a  prediction  of  cloud  motion  12-24  hours  in  advance, 
taking  no  more  than  6  hours  from  initiation  to  completion.] 

A  companion  paper  [3]  will  provide  the  technical  details  of  our  approach  to  uncertainty  quantification  for  this 
problem.  We  provide  here  an  overview.  To  discuss  uncertainty  and  its  effects  in  models,  we  distinguish  between 
parametric  uncertainty  and  stochastic  forcing.  Mathematically,  these  two  kinds  of  uncertainty  are  accounted  for 
differently.  However  there  is  not  a  perfect  demarcation  between  the  way  we  treat  these.  For  example,  we  cannot 
know  the  input  windfield  perfectly  -  at  every  point  in  the  atmosphere  at  every  instant  of  time  and  we  may  choose  to 
model  this  uncertainty  as  a  stochastic  forcing  of  a  given  size  at  each  computational  location  in  space-time.  Following 
arguments  in  [24]  that  probability  theory  is  the  right  setting  to  describe  our  lack  of  information,  we  adopt  here  a 
probabilistic  framework. 
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Instead  of  solving  for  the  ash  concentration  at  a  given  time  and  location  (based  on  integration  of  puf  f  particles) 
we  treat  the  concentration  at  a  given  location  and  time  as  a  random  variable  x^.  The  time  evolution  of  the  vector 
of  x^s,  is  given  by  a  stochastic  differential  equation  (which  should  be  thought  of  as  generalization  of  Eq.  (1)  to  include 
the  physics  of  the  bent -puf  f  coupling): 


X&+1  =  f (tk,  xk,  0)  +  rjk ,  x(t0)  =  x0  (2) 

In  this  equation,  0  represents  uncertain  but  time-invariant  system  parameters  such  as  the  coefficient  of  gravity  fallout, 
vent  diameter  or  vent  velocity.  The  random  function  r\k  represents  stochastic  forcing  terms,  such  as  turbulent  diffusion, 
or  any  stochastic  effect  that  might  be  added  to  the  windfield.  These  random  functions  may  be  modeled  as  a  Gaussian 
process  with  specified  correlation  function  Q^.  The  nominal  initial  position  of  the  center  of  the  particles  is  given 
by  xo,  which  may  also  be  uncertain.  The  total  uncertainty  associated  with  the  state  vector  xk  is  characterized  by  the 
probability  distribution  function  (pdf)  p(tk,  xk,  0)  whose  time  evolution  characterizes  our  uncertainty.  Using  the  full 
pdf  has  many  advantages,  especially  the  ability  in  a  DDDAS  framework  of  integrating  data  and  prediction  easily  using 
Bayesian  approaches. 

Several  approximate  techniques  are  commonly  used  to  approximate  the  state  pdf  evolution  [25,  26],  the  most 
popular  being  Monte  Carlo  (MC)  methods  [27],  Gaussian  closure  [28],  Equivalent  Linearization  [29],  and  Stochastic 
Averaging  [30].  In  addition,  a  Gaussian  Process  approach  to  solve  nonlinear  stochastic  differential  equations  has  been 
proposed  in  Ref.  [32].  All  of  these  algorithms  except  MC  methods  are  similar  in  several  respects,  and  are  suitable 
only  for  linear  or  moderately  nonlinear  systems,  because  the  effect  of  higher  order  terms  can  lead  to  significant  errors. 
Monte  Carlo  methods  require  extensive  computational  resources  and  effort,  and  become  increasingly  infeasible  for 
high-dimensional  dynamic  systems  [33]  and  the  limited  and  time  bound  response  needed  in  DDDAS  settings. 

The  next  two  subsections  discuss  our  planned  methods  for  solving  the  time  evolution  of  state  pdf  for  systems  that 
include  stochastic  forcing  terms  and  parametric  uncertainty. 

3.2.  Parameter  Uncertainty  Propagation 

The  bent-  puff  model  is  subject  to  parametric  uncertainty.  The  propagation  of  uncertainty  due  to  time-invariant 
but  uncertain  input  parameters  can  be  approximated  by  a  generalized  polynomial  chaos  (gPC),  originally  due  to  Xiu 
and  Karniadakis  [34].  Roughly  speaking,  gPC  postulates  a  separation  of  random  variables  from  deterministic  ones 
in  the  solution  algorithm  for  a  differential  equation.  The  random  variables  are  expanded  as  a  polynomial,  those 
polynomials  being  associated  with  the  assumed  probability  distribution  for  the  input  variables  (for  example,  Hermite 
polynomials  for  normally  distributed  parameters,  or  Legendre  for  uniformly  distribution).  Galerkin  collocation  is 
used  to  generate  a  system  of  coupled,  deterministic  differential  equations  for  the  expansion  coefficients.  The  gPC 
collocation  step  fails  when  applied  to  problems  with  non-polynomial  nonlinearities,  and  can  produce  non-physical 
solutions  when  applied  to  hyperbolic  equations.  Non-intrusive  spectral  projection  (NISP)  or  stochastic  collocation 
methods  can  overcome  these  difficulties  [35-37].  We  will  use  a  different  formulation  of  the  NISP  idea  [38]  that  we 
call  polynomial  chaos  quadrature  (PCQ).  PCQ  replaces  the  projection  step  of  NISP  with  numerical  quadrature.  The 
resulting  method  can  be  viewed  as  a  Monte-Carlo-like  evaluation  of  (2),  but  with  sample  points  selected  by  quadrature 
rules. 

The  bent-  puff  model  contains  both  stochastic  forcing  and  parametric  uncertainty.  We  assume  the  paramet¬ 
ric  uncertainty  is  independent  of  the  forcing  terms,  and  propose  wrapping  the  PCQ  quadrature  evaluation  around 
an  AGMM  solution  step.  For  a  fixed  value  of  parameter  0  =  0;,  the  AGMM  provides  an  output  distribution  ac¬ 
counting  for  stochastic  forcing.  Integration  (by  PCQ)  over  the  uncertain  inputs  determines  the  posterior.  Expressed 
mathematically, 


P(*k+ 1) 


J  p(t,xk+i\Q)dp(Q)  ~  y^toqp(t ,xk+\\®q)p(®q) 
J  q 

w[+lN(xk+ 1  \n‘k+vP *+1,09)j />(©,) 


(3) 


In  this  expression,  cjq  are  the  PCQ  weights  assigned  to  the  quadrature  points  Gq. 
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3.3.  Uncertainty  Propagation  due  to  Stochastic  Forcing 

Consider  first  Eq.  (2)  for  fixed  system  parameters  ©o,  and  a  stochastic  forcing  77.  In  the  context  of  the  bent-  puff 
model,  this  assumption  is  akin  to  fixing  all  input  parameters  but  allowing  for  a  random  component  the  windfield. 
The  time  evolution  of  the  pdf  from  time-step  k  to  time-step  k  +  1  is  given  by  the  Chapman- Kolmogorov  equation 
(CKE)  [39-41]: 

p(tk+ 1,  x*+i|0o)  =  J  p(tk+i,xk+1\tk,  xk,  &o)p(tk,  Xfc|0o)dxfc  (4) 

Here  pfa+i,  x^+i|4,  x^,  ©0)  is  the  state  transition  corresponding  to  the  pdf  for  the  process  noise  77^  [42].  The  CKE  is 
a  formidable  equation  to  solve  because:  (/)  the  pdf  must  remain  positive,  (ii)  the  pdf  must  maintain  unit  volume,  and 
(iii)  the  domain  of  p  changes  in  time. 

The  integral  in  Eq.  (4)  is  over  the  state  space  x^,  so  it  is  high  dimensional.  Traditional  numerical  approaches  which 
discretize  the  space  in  which  the  pdf  lies  suffer  from  the  “curse  of  dimensionality”.  To  overcome  this  obstacle,  we 
will  use  a  recently  developed  adaptive  Gaussian  mixture  model  (AGMM)  in  conjunction  with  information  theoretic 
measures  [43,  44]  to  accurately  solve  the  CKE.  The  key  idea  of  the  AGMM  is  to  approximate  the  state  pdf  by  a 
finite  sum  of  Gaussian  density  functions  whose  mean  and  covariance  are  propagated  from  one  time-step  to  the  next 
using  linear  theory.  The  weights  of  the  Gaussian  kernels  are  updated  at  every  time-step,  by  requiring  the  sum  to 
satisfy  the  CKE  [45].  When  properly  formulated,  the  mixture  problem  can  be  solved  efficiently  and  accurately  using 
convex  optimization  solvers,  even  if  the  mixture  model  includes  many  terms.  This  methodology  effectively  decouples 
a  large  uncertainty  propagation  problem  into  many  small  problems.  As  a  consequence,  the  solution  algorithm  can  be 
parallelized  on  most  high  performance  computing  systems.  Sec. 4  below  briefly  descibes  preliminary  results,  where 
an  ensemble  of  bent-  puff  computations  were  run  in  parallel.  Finally,  a  Bayesian  framework  can  be  used  on  the 
AGMM  structure  to  assimilate  (noisy)  observational  data  with  model  forecasts  [46,  47]. 

3.4.  Data  Assimilation  for  Reducing  Uncertainty 

Of  course  using  any  sensor  data  that  might  become  available  to  correct  and  refine  the  dynamical  model  forecast 
will  reduce  the  uncertainty  of  predictions.  This  model-data  fusion  process  is  well  documented  in  meteorology  appli¬ 
cations  where  a  variety  of  data  assimilation  methods  are  used  operationally  to  improve  the  quality  of  the  Numerical 
Weather  Prediction  (NWP)  forecast[48-51].  Several  recent  studies  additionally  used  conventional  filtering  approaches 
to  incorporate  observations  of  a  chemical  species  into  a  dispersion  model[52,  53]  and  discuss  how  the  nonlinearities 
in  T&D  models  and  sparse  measurement  data  cause  problems  in  capturing  non-Gaussian  state  pdf.  In  general,  find¬ 
ing  an  optimal  filter  is  an  infinite  dimensional  problem  [54,  55].  Finite  dimensional  filters  have  been  derived  for  a 
certain  nonlinear  problems  [56-59].  For  more  general  but  low-order  nonlinear  systems,  particle  filter  methods  hold 
promise  [60]. 

Standard  Bayesian  tools  will  be  employed  to  incorporate  observational  data.  Given  a  prediction  of  the  state 
variable  x^,  standard  Bayesian  algorithms  assume  a  sensor  model  h  to  obtain  the  measurement,  y^  =  h(4,x^)  +  v^. 
Here  v^  is  the  measurement  noise  with  prescribed  likelihood  function  p( y^|x^).  Using  the  dynamic  state  evolution 
sketched  above  as  a  forecasting  tool,  the  state  pdf  can  be  updated  using  the  Bayes’  rule  and  measurement  data: 

IV  s  p(yk\xk)p(xk\Yk-i) 

P  f  p(y*|Xfc)p(x*|Y*_i)dx* 

Here,  p(x^|Y^_i)  represents  the  prior  pdf  (the  computed  solution  of  the  CKE),  p( y^x^)  is  the  likelihood  that  we  observe 
y k  given  the  state  x^  and  p(x^|Y^)  represents  the  posterior  pdf  of  x&.  For  a  fixed  set  of  data  and  underlying  probability 
model,  Maximum  a-Posteriori  (MAP)  picks  the  values  of  the  state  estimates  that  make  both  the  data  and  model  “more 
likely”  than  any  other  values  of  the  state  estimates  would  make  them.  This  may  be  of  considerable  use  for  decision 
support  in  situations  in  which  multiple  strong  and  roughly  equal  modes  appear  in  the  posterior  estimate.  Neither  the 
conditional  mean  nor  the  maximum  likelihood  estimate  as  the  sole  reported  hypotheses  are  sufficient  to  formulate 
effective  decisions  in  this  case  [61]. 
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3.5.  Model  Hierarchy 

The  use  of  the  modestly  complicated  puf  f  code  and  the  availability  of  other  models  of  differing  fidelities  opens 
up  the  possibility  of  using  multiple  models  with  the  higher  resolution  model  acting  as  an  effective  alternate  provider 
for  gaps  in  the  sparse  data  from  observations.  We  will  focus  these  runs  on  well  defined  parameter  ranges  where  the 
basic  puff  code’s  output  is  deemed  to  be  poor  based  on  observations.  Our  preliminary  work  (see  below  4)  and  recent 
literature [62]  indicates  that  the  puff  code  is  sensitive  to  the  choice  of  number  of  particles.  However,  the  cost  of  the 
computation  with  very  large  numbers  of  particles  (107  as  opposed  to  the  more  common  105  )  is  very  high.  We  will 
selectively  use  a  few  of  these  high  particle  runs  to  provide  a  high  fidelity  source  where  we  notice  large  discrepancy 
between  the  satellite  and  bent -puf  f  outputs. 

3.6.  Observation  and  Testing 

We  are  fortunate  to  have  access  to  a  wealth  of  data  from  the  2010  Eyjafjallajokull  eruption  and  the  subsequent 
cloud  dispersion.  The  daily  windfields  for  the  important  period  between  April  14-20  are  known,  satellite  images 
tracking  the  cloud  are  available,  as  is  plume  data  (for  example,  there  is  Meteosat  SEVIRI  data  on  ash  loading  and 
plume  height,  and  SEVIRI  and  CALIPSO  data  showing  the  motion  of  the  cloud).  This  event  will  provide  a  rigorous 
test  of  the  proposed  uncertainty  methodology,  and  in  particular  the  sensitivity  of  simulation  results  to  differing  plume 
and  windfield  inputs. 

Webley  et  al.  [63]  suggest  the  possibility  of  running  VATD  models  using  retrieved  ash  cloud  parameters  from 
satellite  remote  sensing  data.  We  will  assess  this  potential.  We  will  also  generate  ‘synthetic’  satellite  imagery  from 
the  puff  model  simulations,  using  the  ash  concentrations,  optical  depths  of  the  clouds  and  particle  sizes  to  generate 
equivalent  reverse  absorptions  signals  and  opaque  cloud  temperatures.  This  process  is,  essentially,  an  inverse  of  the 
volcanic  ash  retrieval  process  developed  by  Wen  and  Rose  [64]  and  Pavolonis  et  al.  [65].  Here,  we  will  use  the 
VATD  concentration  maps  to  determine  the  effect  the  ash  would  have  on  the  temperature  of  the  cloud  as  seen  by  the 
satellite.  From  this  we  can  compare  datasets  and  examine  the  satellite  retrieval  limits  and  reverse  absorption  detection 
capabilities.  This  synthetic  imagery  from  the  puff  model  will  then  allow  statistical  comparisons  between  the  different 
datasets,  to  ascertain  which  parameters  have  the  greatest  effect  on  the  model  simulation  uncertainties. 


4.  Preliminary  Results 

A  first  step  in  this  research  effort  is  reported  in  [66],  where  four  bent  inputs  (vent  radius,  vent  velocity,  mean  grain 
size,  and  grain  size  variance)  were  considered  as  uncertain,  bent  outputs  fed  into  an  ensemble  of  puff  simulations, 
and  PCQ  was  employed  as  a  solution  methodology  (see  additional  details  in  companion  paper  [?  ]. 

Following  runs  of  bent  at  Clenshaw-Curtis  quadrature  points,  each  bent  output  is  then  propagated  through  puff, 
which  was  then  run  for  a  real-time  period  of  five  days.  The  outputs  from  puff  were  then  combined  to  produce  the 
ensemble  by  applying  the  appropriate  weight  to  each  deterministic  bent-puff  run. 

Numerical  methodology  of  the  type  used  here  must  always  be  tested  for  consistency,  i.e.,  independence  from 
discretization  parameters  or  sample  size.  As  we  seek  to  develop  an  ensembling  method  that  could  be  implemented 
in  a  finite  computing  time,  we  also  need  to  minimize  the  number  of  parcels  and  runs.  For  this  set  of  simulations 
we  compared  the  simulation  outputs  for  94  and  134  samples  (quadrature  points)  and  the  results  indicate  that  94  runs 
yields  substantially  the  same  results  as  134  runs.  The  comparison  of  using  105,  4  x  106  and  107  particles  in  the  puff 
simulation  also  indicated  that  the  choice  of  4  x  106  particles  was  adequate  for  our  purposes,  and  consistent  with  the 
findings  of  others  [62] . 


5.  Conclusions  and  Future  Directions 

Model  output  was  compared  with  Meteosat-9  retrievals  of  top  ash  plume  height.  Volcanic  ash  was  identified  in 
the  satellite  data  using  the  methodology  described  in  [67]  and  [68].  The  ash  cloud  height  were  retrieved  using  an 
optimal  estimation  approach  [69]  and  [70].  All  microphysical  assumptions  used  in  the  retrieval  are  described  in  [71]. 
In  general,  plume  location  is  quite  close  to  that  shown  in  the  data,  with  little  variance.  There  are  greater  divergences 
between  model  and  data  in  the  eastern  portion  of  the  plume  where  Meteosat-9  data  estimate  a  much  smaller  plume 
area  than  does  the  model.  Estimating  probabilities  by  computing  the  coefficients  in  the  polynomial  approximation  of 
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(c)  Ash  Top  Height  (16r/7  April,  1200hrs) 


(d)  Ash  Top  Height  (16th  April,  1800hrs) 


(e)  Probability  of  ash  based  on  PCQ  based  ensemble  at  \6th  April,  1800hrs 

Figure  2:  Meteosat-9  SEYIRI  ash  top  height  data  compared  with  simulation  output  for  a  94  run  implementation  of  PCQ  sampling 
of  the  input  space.  106  ash  particles  were  used  in  the  simulation.  Probability  of  ash  based  on  PCQ  based  ensemble  at  \6th  April, 
1800hrs  from  [3];  Inner  contour  Prob  >  0.7,  outer  contour  Prob  >  0.2 
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the  pdf  from  the  PCQ  methodology  indicates  that  all  observed  ash  is  within  the  0.2  probability  contours  and  most  is 
within  the  0.7  contour.  It  appears  from  the  literature  that  this  is  the  first  time  such  quantitative  information  has  been 
available  from  simulation  albeit  in  hind  cast  mode.  In  short,  preliminary  results  shows  that  PCQ  scheme  yields  a 
valid  estimate  of  mean  and  standard  deviation  of  ash  parameters  for  a  relatively  small  sample  size  of  order  [94].  The 
question  arises  however,  whether  even  PCQ  can  provide  probabilistic  estimates  useful  in  Volcanic  Ash  Advisories  in 
a  reasonable  time.  This  question  clearly  needs  to  be  investigated  further. 
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